/***************************************************************************
 *   Copyright (C) 2008 by Ivan Blagoev Topolsky   *
 *   ivan.topolsky@genebio.com   *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU Library General Public License as       *
 *   published by the Free Software Foundation; either version 3 of the    *
 *   License, or (at your option) any later version.                       *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU Library General Public     *
 *   License along with this program; if not, write to the                 *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <iostream>
#include <iomanip>

// BOOST
#include <boost/regex.hpp>
#include <boost/lexical_cast.hpp>

// Captain-bol
#include <spectrumpeakrefdualptr.h>
#include <chargemapper.h>

// OLE/COM support
#include "comutil.h"
#include "disphelper.h"

#include "filterparse.h"
#include "arthur.h"

namespace xeno {

namespace merlin {

using namespace std;	// STL
using namespace bol;	// Captain bol
	
// FIXME warning this isn't bytesex-portable - on the other hand, I don't think that we will ever run this code on anything but x86-compatible chips
// although, for the record, the "Darwine" project enables one to use Windows code through Wine/winelibs
// with the help of an x86 emulation layer by special "*nix-inside-*nix" mode of QEMU 
// (= doesn't emulate a full x86 machine, only helps calling x86 code in *nix .so libraries from another *nix. In in this case, the .so is Winelib).
// ... SO IN FACT THE WHOLE THING COULD EVEN BE PORTED TO MAC OS X ON PPC !!! ... [Yann would be happy ;-) ]
typedef union {
	long long refid;
	struct {
		unsigned short scanid;
		unsigned short detectid;
	} th;
} REFID;

const string arthur::default_runname_template ( DEFAULT_RUNNAME_TEMPLATE );

/** 
-=- !!!XDK!!! -=-

Explanations about the XDK and how to interface Xcalibur will be interwoven in the code's comments.

General overview :

raw file structure :		method using :	comments about method :
XRaw				initfile	loads file, initiate the main OLE/COM object
|--XRawInfo			 :
|--XSample			samplelsit	lists all general information about a sample
|--XDetectors			dosample	get all MS data about the sample by calling the next methods
|  +--XDetector			 :	(also used by dospectra_all to access the filters)
|--XFilters			 :
|  +--XFilter		class	filterparse 	all information about scan and spectra is in a simple string
+--XSpectra			dospectra [_all]	get all the spectra
   |					dospectra - called by dosample, once for each different filter
   |					dospectra_all - called just once and process all filters
   +--XSpectrum			do_single_spectrum_ms#	get the data of 1 sample
      |--XParentScans		 :
      |  +--XParentScan		 :		data is cached an used by 'buildptr' to make link between spectra
      +--Win32::SafeArray	getdata		the peak raw data is stored in a standard OLE/COM SafeArray

NOTE to access OLE/COM we use a slightly modified version of "disphelper" - http://disphelper.sourceforge.net/
- it uses late binding, calling methods by text discribing their name 
  (no DLL needed at compile time, can be compiled by Linux)
- it has a printf-like syntax :
	dhCheck( dhCreateObject(L"MSXML.DOMDocument", NULL, &xmlDoc) );
	dhCheck( dhPutValue(xmlDoc, L".Async = %b", FALSE) );
	dhCheck( dhCallMethod(xmlDoc, L".Load(%s)", szURL) );
	dhGetValue(L"%s", &RSSItem.szTitle,       xmlNode, L".selectSingleNode(%S).text", L"title");
- my modification includes supporting "arguments as references" using a syntax like "%&d"
  (the original disphelper only supported references packaged by the user in VARIANT objects)

*/


/*
 * ENUMs from XDK
 */
 

/// Whenever there's an enum, XDKNAMES maps a nice name for print
XDKNAMES::XDKNAMES() {
	// following code was programmatically generated :
	XDetectors[XNo_Device] =  "No Device";
	XDetectors[XMS_Device] =  "MS Device";
	XDetectors[XMS_Analog_Device] =   "MS Analog Device";
	XDetectors[XAnalog_Device] =      "Analog Device";
	XDetectors[XPDA_Device] =         "PDA Device";
	XDetectors[XUV_Device] =  "UV Device";

	SampleTypes[XSampleUnknown] =      "Unknown";
	SampleTypes[XSampleBlank] =        "Blank";
	SampleTypes[XSampleQC] =   "QC";
	SampleTypes[XSampleStdClear] =     "StdClear";
	SampleTypes[XSampleStdUpdate] =    "StdUpdate";
	SampleTypes[XSampleStdBracket] =   "StdBracket";
	SampleTypes[XSampleStdBracketStart] =      "StdBracketStart";
	SampleTypes[XSampleStdBracketEnd] =        "StdBracketEnd";
	SampleTypes[XSampleProgram] =      "Program";

	/*
	[XBarcodeNoRead] =      "NoRead";
	[XBarcodeRead] =        "Read";
	[XBarcodeUnreadable] =  "Unreadable";
	[XBarcodeError] =       "Error";
	*/
}



/**
* -=- !!!XDK!!! -=-
* One peculiarity of XCalibur is that it doesn't describe samples with a simple textfield (name) 
* But instead has tons of properties describing it,
* including 5 fields (User 1 to 5) whose meaning a left free to the user to decide
*
* The crazy amount of sample properties follows :
*/
const wchar_t *arthur::fields_str[] = {	// String properties. Return BSTR
	L"SampleId",	L"SampleName",	L"Comment",	L"Barcode",	L"Level",	L"Position",
	/*L"Path",	L"FileName",	L"InstMeth",	L"ProcMeth",	L"CalFile",*/	NULL
},
*arthur::fields_dbl[] = {	// Return double
	L"DilFactor",	L"SampleWt",	L"SampleVol",	L"InjVol",	L"ISTDAmt",	NULL
},
*arthur::fields_num[]= {	// return some integers
	/* long: */ L"Revision", /* short : */ L"RowNumber", NULL
	/*Enums :	XSampleTypes SampleType,	XBarcodeStatus BarcodeStatus */
},
*arthur::fields_enu[]= {	// return some some ENUMs
	L"SampleType",	/*XBarcodeStatus BarcodeStatus, */ NULL
};



/**
 * Creates an nice string name out of a template
 *  We will use string substitution de replace "#Fields" with actual value.
 *
 * @param name	the (template) name that we are building
 * @param label the name of the Field (without the heading '#')
 * @param val 	the value to substitute
 */
inline string& arthur::nfo(string& name, const string& label, const string& val) {
	using namespace boost;
	using namespace boost::regex_constants;

	const regex re("#" + label, literal | icase);
	return name = regex_replace(name, re, val, match_default);
}

inline string& arthur::nfo(string& name, const wchar_t *label, const string&  val) {
	string needle;
	wcstostring(needle, label);

	return nfo(name, needle, val);
}



/**
 * This template is used to process all Sample properties which return a given TYPE
 * - It prints them nicely on the console
 * - If a template string is provided, it will replace all '#fields' using the bove 'NFO' method
 *
 * typename type : the type
 * @param name	non-NULL: string template. NULL: print only
 * @param fmt	scanf/printf-like format used by disphelper (%d, %e, %s, ...)
 * @param sample	XSample object whose properties we are printing
 * @param f	point to one of the 'fiedls_???' list of property names above
 * @param mptr	non-NULL : one of the above map for enums which converts the raw propertie values to nice text names, NULL: put the value as-is 'name'
 */
template <typename TYPE> inline void arthur::shownfo(string *name, const wchar_t *fmt, IDispatch *sample, const wchar_t **f, map<TYPE,string>*mptr) {
	TYPE data;
	string val;
	for(;*f;f++) {
		dhCheck( dhGetValue(fmt, &data, sample, *f) );
	
		val =  (mptr)
			? (*mptr)[data]	// print the data as is
			: (boost::lexical_cast <string> (data)); // use a map-pointer to put a nice string on screen
	
	#if defined(__GNUC__) && (__GNUC__<4)
		// MinGW is based on an older version of GCC which picks up the wrong operator if not forced :
		::operator<< ((ostream&)(cout << " - "), (*f)) << " :\t" << val << endl;
	#else
		cout << " - " <<  (*f) << " :\t" << val << endl;
	#endif	
		if (name) 
			nfo(*name, *f, val);
	}
};




/**
 * this function takes care of the first steps of handling raw files :
 * - getting the necessary OLE/COM object
 * - setting up a RunLcmsms object to hold the results
 * 
 * @param ifname Filename to use
 * @param ptr Write to a specific RunLcmsms instance. Other wise, if NULL, we allocate one.
 */
void arthur::initfile(const char *ifname, RunLcmsms *ptr) {
	// reset the tables
	// as the spectra are all represented using USHORT in XCalibur, I don't except to have more than 64k spectra to handle
	parentmap.clear(); parentmap.resize(USHRT_MAX,0);
	parentsp1.clear(); parentsp1.resize(USHRT_MAX,NULL);
	parentsp2.clear(); parentsp2.resize(USHRT_MAX,NULL);
	extractmap.clear(); extractmap.resize(USHRT_MAX,false);
	
	// *** EXPORT CAPTAINBOL
	
	// allocate an instance if needed...
	run = ptr ? ptr : new RunLcmsms();

	//
	// Open the file
	//

	/** 
	 * -=- !!!XDK!!! -=-
	 * XDK presents an API for accessing RAW files which consist of a RAW-file reader (XRaw)
	 * which presents the data in a structured content.
	 * See "Raw File Hierarchy" in the XDK documentation :
	 *	C:/Xcalibur/help/xdkhelp/rawfile_hierarchy.htm
	 * (unlike, for example, Analyst Cookbook which present an API more oriented toward piloting Analyst
	 * and presents data in a way which closely maps to various phase of the pipeline)
	 *
	 * in Xcalibur, each RAW file contains one and only one single sample, and all data pertaining it
	 * the sample is the first-grade citizen and everything is structured around it.
	 *
	 * Note about classes :
	 * - every count is 1-based and everything is counted consistently using 1...n like in Basic (not 0...n-1 like in C)
	 * - the classes aren't grouped in a common name space like "Xcalibur.{class}.1" or "XDK.{class}.1"
	 *  instead the classes use the same name for name space "{class}.{class}.1"
	 * - several classes have "~Read" equivalent which are read only.
	 */
	//DISPATCH_OBJ(raw);
	try {	dhCheck( dhCreateObject(L"XRaw.XRaw.1", NULL, &raw) );
	} catch (string errstr) {
		cerr << "Fatal error - can't create the main XCalibur object (XRaw)" << endl
			<< endl << "Is the library correctly registered ?" << endl
			<< "- either run Xcalibur's installer" << endl
			<< "- or manually run 'regsvr32 {drive}:\\{path}\\XcaliburFiles.dll'" << endl
			<< "(as explained in the 'Xcalibur Files Frequently Asked Questions (FAQ)' section of the XDK)" << endl
			<< endl << "details:" << endl << errstr << endl;
		exit(2);
	}

	cout << "opening " << ifname << " ... ";
	try {	dhCheck( dhCallMethod(raw, L"Open %s", ifname) );
	} catch (string errstr) {
		cerr << "Fatal error - can't open the file" << endl
			<< "details:" << endl << errstr << endl;
		exit(2);
	}	
	cout << "done" << endl;

	//
	// Get information 
	//
	
	/**
	 * -=- !!!XDK!!! -=-
	 * XRawInfo info holds generic info about file :
	 * - Number & types of detectors
	 * - Title of each user field
	 *
	 * XSample holds information about the sample analyzed in this file :
	 * - a dozen of fields specifying data
	 * - 5 user defined fields (each field described in XRawInfo - see above)
	 */
	//DISPATCH_OBJ(rawnfo);
	/*int*/ numdetect = 0;	// number of detectors
	
	try {
		dhCheck( dhGetValue(L"%o", &rawnfo, raw, L"RawInfo") );
		
		//
		// Get number of detectors
		//
		dhCheck( dhGetValue(L"%d", &numdetect, rawnfo, L"DetectorCount") );
	}  catch (string errstr) {
		cerr << "Fatal error getting info - details:" << endl << errstr << endl; 
		exit(2);
	}
	
	
	// Normally, run flow continues with the next function
	// they are currently split to allow separate calling 
}

/**
 * normally this step comes right after the first
 * - it will display all information about a sample
 * - given a template, it will replace all "#Fields" with proper values
 * 
 * @param nametmpl Template out of which to build the name
 *
 * Also relies on the XRawInfo stored in 'rawnfo'
 */
void arthur::samplelist(const string* nametmpl) {
	//
	// Show all information about sample
	//
	try {
		WITH(sample, raw, L"Sample") {
			string desc (nametmpl ? *nametmpl : ""); // holds the description of the run
				// "#Fields" will be replaced with corresponding value

			cout << "Sample:" << endl;

			// Ouput standard information about sample
			shownfo<XSampleTypes>	((nametmpl?&desc:NULL), L"%d", sample, fields_enu, &XDKNames.SampleTypes);
			shownfo<char *>	((nametmpl?&desc:NULL), L"%s", sample, fields_str); // FIXME there's a memory leak (string should be freed)
			shownfo<double>	((nametmpl?&desc:NULL), L"%e", sample, fields_dbl);
			shownfo<int>	((nametmpl?&desc:NULL), L"%d", sample, fields_num);

			// Ouput user-defined fields about sample
			for (int i = 1; i <= 5; i++) {
				char *name = NULL, *data = NULL;
				dhCheck( dhGetValue(L"%s", &name, rawnfo, L"UserLabel %d", i) ); // FIXME memory leak happens here...
				dhCheck( dhGetValue(L"%s", &data, sample, L"UserText %d", i) );  // ...if exception is thrown there
				cout	<< "[" << i << "] " << name << " :\t"
					<< data << endl;
				
				// if we got a template, replace the "#fields"
				using namespace boost;
				if (nametmpl)	nfo(desc,
					"User" + lexical_cast <string> (i),
					lexical_cast <string> (data));
				
				SAFE_FREE_STRING(data);
				SAFE_FREE_STRING(name);
			}
			
			// if we have a template store the name from the template
			if (nametmpl)
				// *** EXPORT CAPTAINBOL
				run->setDescription(desc);
		} END_WITH_THROW(sample);
	} catch (string errstr) {
		cerr	<< "Error getting info - details:" << endl 
			<< errstr << endl
			<< "Trying to continue anyway with partial data..." << endl;
	}

	// Normally, run flow continues with the next function
	// they are currently split to allow separate calling 
}


/**
 * process all MS data pertaining the current file's sample
 * normally this step comes right after the above one
 *
 * @param allatonce false : load filters one by one	true : loads all sample at once (might be fast in Wine)
 * @param skipms1	if true the peak data for MS precusor spectra will be skipped and only structure will be read
 * @param skipms2		idem for MS^2 fragment spectra
 * @param dropms1	if true completely skip everything about MS precursors, only process MS^2 fragment
 * @param dropms2		idem for MS^2 fragment spectra
 * @param list 	if non-NULL: a list of 64bits refID's whose spectra will be processed any way, regardless of the skipms# boolean
 * The trick is to usually set skipms# to true, and only load the needed data using 'list'
 */
void arthur::dosample(const bool allatonce, const bool skipms1, const bool skipms2, const bool dropms1, const bool dropms2, ExtractList* reflist) {
	if (numdetect == 0) {
		cerr << "Run initfile first" << endl;
	}


	//
	// Proceed Detectors
	//

	/**
	 * -=- !!!XDK!!! -=-
	 * a XDetector represents measures that were done on the sample using some specific instrument
	 * normally, we're going to encounter only 1 detector : some mass spectrometer
	 *
	 * but, a XRaw can hold several different measures done on the same sample
	 * the classical example is first using a UV light source to read the output of a Liquid Chromatography
	 * and then analyzing with a mass spectrometer
	 *
	 * there are even provitions in the .RAW format to have several detector of the same type
	 * theoretically, one could imagine using a tee at the exit of the Chromatography to split the elute
	 * - part of the elute goes into an ESI and is analyzed online
	 * - part of the elute goes to a plotter-robot which plots on a MALDI plate for further offline analysis
	 * all measure would be stored as separate XDetector in the XRaw
	 * in practice this kind of setup is seldom used (because of the costs)
	 */
	cout << numdetect << ((numdetect > 1) ? " detectors" : " detector") << endl;
	for (int d = 1; d <= numdetect; d++) {
		initextract_map(d, reflist);

		
		/* int */ numfilter = 0;
		// DISPATCH_OBJ(detect);
		try {
			bool isaMS = false;
			XDetectorTypes type;
			
			//
			// Detector type
			//
			dhCheck( dhGetValue(L"%d", &type, rawnfo, L"DetectorType %d", d) );
			cout << "Detector " << d << " is type " << type << " ";
			
			switch (type) {
				case XNo_Device:	cout << "No Device";	break;
				case XUV_Device:	cout << "UV Device";	break;
				case XPDA_Device:	cout << "PDA Device";	break;
				case XAnalog_Device:	cout << "Analog Device";	break;
				case XMS_Analog_Device:
					cout << "Analog "; // chain..
				case XMS_Device:
					cout << "MS Device";
					dhCheck( dhGetValue(L"%o", &detect, raw, L"Detector %d,%d", type, d) );
					isaMS = true;
					break;
			}
			
			if (! isaMS) {
				cout << "...skipping" << endl;
				continue;
			}
		} catch (string errstr) {
			cerr	<< "Error getting info detector" << endl 
				<< errstr << endl
				<< "Trying to continue with next detector..." << endl;
			continue; // jump straight to next detector. Do not attempt to get filters
		}


		/**
		 * -=- !!!XDK!!! -=-
		 * Other data :
		 * One day some application may be interested in :
		 * - XInstrument which hold data about instrument
		 * - XTune which hold data about the instrument tuning
		 */


		//
		// Filters
		//
		
		/**
		 * -=- !!!XDK!!! -=-
		 * Inside an XDetector all measurments (the spectra) 
		 * aren't just all sequentially available ( NOT Spectrum[num] )
		 * instead, spectra that share the same detector settings (for an MS: selected M/Z, collision energy, etc) 
		 * are grouped together by filters. (i.e.: Spectrum[filter][num] )
		 *
		 * helper functions can find inside which filter a given spectrum is.
		 *
		 * also, by using "0" as the filter number when manipulating spectra collection, 
		 * one gets to use the whole spectra at once
		 * (which usually fucks up the 32bits memory space barrier usef by Windows applications and OLE/COM)
		 *
		 * an XFilter hold infromation about the detector setting
		 * SADDLY, all that data is only available as a string
		 * (though there's a method available to check validity of a string)
		 */

		bool filterdone = false;
		if (! allatonce) try { 
			WITH(filtercol, detect, L"Filters") {
				/*int*/ numfilter = 0;
				dhCheck( dhGetValue(L"%d", &numfilter, filtercol, L"Count") );
				cout << " and has " << numfilter << " filters" << endl;
				
				
				/**
				 * -=- !!!XDK!!! -=-
				 * The data about filters is stored inside a text string
				 * We first use a parser to process it
				 */
				filters.resize(numfilter+1, NULL);
				
				for (int f = 1; f <= numfilter; f++) {
					try {
						WITH1(filter, filtercol, L"Item %d", f) {
							char *str = NULL;
							dhCheck( dhGetValue(L"%s", &str, filter, L"Text") );
							cout << "\tfilter " << f << " " << str;
							
							// store the processed string in the vector which store together all the filters.
							filters[f] = new FilterParse (str);
							
							cout << endl;
							SAFE_FREE_STRING(str);
						} END_WITH_THROW(filter);
					} catch (string errstr) {
						cerr	<< "Error getting filters collection" << endl 
							<< errstr << endl
							<< "Trying to continue with next filter..." << endl;
						continue;
					}
				}
				cout << endl;
				filterdone = true;
				
				/**
				 * Extact spectra based on type of scan
				 */
				for (int f = 1; f <= numfilter ; f++) {
					FilterParse*  curfilt = filters[f];
					
					// Non decoded : not worthy trying
					if (curfilt == NULL) 
						continue;

					// Thankfull FilterParser has now the data in a usable form
					cout << "extracting filter " << f
						<< "\tType " << curfilt->ScanType() << " MS" << curfilt->MSLevel();

					// Process the data according the type...
					if (curfilt->ScanType() == "") {
						// well, obviously *IF* we have a type
						cerr << "No Scan Type could be parsed from string. Skipping !!!" << endl;
						continue;
					} else if (curfilt->ScanType() == "SRM" || curfilt->ScanType() == "CRM") {
						// FIXME build precursor scans out of Reaction monitor
						cout << "\tunsupported, skipping" << endl;
						continue;
					} else if (curfilt->MSLevel() == 2) {
						if (! dropms2) {
							dospectra(&arthur::do_single_spectrum_ms2, d, f, skipms2);
							cout << endl;
						} else
							cout << "\tuser skipping all" << endl;
						continue;
					} else if (curfilt->MSLevel() == 1) {
						if (! dropms1) {
							dospectra(&arthur::do_single_spectrum_ms1, d, f, skipms1);
							cout << endl;
						} else
							cout << "\tuser skipping all" << endl;
						continue;
					}

					// TODO implement other scans
					cout << "\tunknown, skipping" << endl;
				}
			} END_WITH_THROW(filtercol);
		} catch (string errstr) {
			cerr	<< "Error getting filters collection" << endl 
				<< errstr << endl
				<< "Trying to continue with next detector..." << endl;
		}
		
		/**
		 * Alternate code path :
		 * Instead of using filters, we attempt to load everything at once
		 */
		if (allatonce || (filterdone == false))
			dospectra_all(d, skipms1, skipms2, dropms1, dropms2);

		//SAFE_RELEASE(detect);
	}
}



/**
 * This step is called to process all spectra in a given filter
 *
 * Parameter defining Action :
 * @param effector decide wether we'll be processing parent MS or fragment MS^2
 * its a pointer to the corresponding function whitch will take care of the spectrum
 *
 * General parameters :
 * @param d	detector we are currently to processing 
 * @param f	filter we are currently to processing 
 *
 * Parameters how to process the data
 * @param docentroid	FIXME unimplemented - seems to be lacking in object internally of "XSpectrumRead" coclass.
 * @param skip	do we process the data ? or do we skip it
 */
void arthur::dospectra(const EFFECTOR_FUNC_PTR effector, const int d, const int f, const bool skip) {
	FilterParse* curfilt = (f>=0) ? filters[f] : NULL;
	
	if (curfilt) 
		cout << "\tParent " << curfilt->Parent();

	/**
	 * -=- !!!XDK!!! -=-
	 * The "spectra" collection holds all the spectra, obviously.
	 *
	 * Although the spectra are accessed by number (1 to n)
	 * - Spectra have an absolute number. 
	 * - Helper function convert between index and absolute number 
	 * - Helper function convert between index and retention time
	 * (Retention time is in minute inside Xcalibur)
	 *
	 * Those helper functions will help us getting more information
	 * Also we can use them when searching for some specific Scans
	 */
	try {
		/**
		 * -=- !!!XDK!!! -=-
		 * According to FAQ, ***ACCESSING "XSPECTRA" IS SLOW*** 
		 * because the whole collection is calculated and cached
		 *
		 * This is specially obvious under linux (with native elfs built against winelib,
		 * and to a lesser exten when running windows exes with wine).
		 *
		 * thus : 
		 * 1. Cache the "XSpectra" collection  (avoid stacking calls "Detect.Spectra(x).Item(j)")
		 * 2. Don't display collection size if spectra are skiped
		 *
		 * See "dospectra_all" for an alternative using a single call to Spectra
		 */
		WITH1(spectcol, detect, L"Spectra %d", f) { // <- this takes time
			int numspect = 0;
			dhCheck( dhGetValue(L"%d", &numspect, spectcol, L"Count") );
			cout << "\t" << numspect << " spectra" << endl;

			for (int s = 1; s <= numspect; s++) {
				try {
					unsigned short	realidx = 0, // plain simple zero-based index as in C (instead of 1-based as in Basic)
							scannum = 0; // absolute number of the scan (important when searchback a specific number)
					double retentiontime = -1;
					bool docentroid = false;

					cout << "\t\tspectrum " << s;

					/*
 					 * Get information from collection Index (Scan number, retention time, etc.)
					 */
					if (SUCCEEDED( dhCallMethod(spectcol, L"IndexToScanNumber %d,%&hd,%&hd", s, &realidx, &scannum) )) {
						cout << " " << scannum;
					} else {
						cerr << "Warning, can't get an absolute scan number !!!" << endl;
					}
					
					if (SUCCEEDED( dhCallMethod(spectcol, L"IndexToRetentionTime %d,%&hd,%&le", s, &realidx, &retentiontime) )) {
						// Xcalibur gives us minutes, we use seconds 
						retentiontime *= 60.;
						cout << " at " << setprecision(4) << retentiontime << "s";
					} else {
						cerr << "Warning, can't get a retentiontime !!!" << endl;
					}
					
					
					WITH1(spect, spectcol, L"Item %d", s) {
						(this->*effector) (
							d, f, 
							spect, scannum, retentiontime, 
							skip, docentroid
						);
					} END_WITH_THROW(spect);
				} catch (string errstr) {
					cerr	<< "Error getting spectrum collection" << endl 
						<< errstr << endl
						<< "Trying to continue with next spectrum..." << endl;
					continue;
				}
			}
		} END_WITH_THROW(spectcol);
	} catch (string errstr) {
		cerr	<< "Error getting spectra collection" << endl 
			<< errstr << endl
			<< "Trying to continue with next filter..." << endl;
	}
}






/**
 * Alternative code-path
 * this alternatives attempts to load all spectra at once
 *
 * this might be slightly faster under Linux as Wine's support for OLE/COM has some overhead.
 * perhaps some future Wine version will be better performing
 *
 * this is *REALLY* hard to achieve in 32bits Windows, because Physical Adress Extension (PAE)
 * is disabled on anything but the most expensive Server Pro Deluxe license
 * limiting the OS to a 4GB memory capacity, and restricting the user address space to 2G
 * (the rest of the address space is allocated to the kernel)
 *
 * And given the fabulous support for 64bits Windows ...
 *
 * General parameters :
 * @param d	detector we are currently to processing 
 * @param f	filter we are currently to processing 
 *
 * Parameters how to process the data
 * @param skipms#	do we process the data ? or do we skip it
 */
void arthur::dospectra_all(const int d, const bool skipms1, const bool skipms2, const bool dropms1, const bool dropms2) {
	/**
	 * -=- !!!XDK!!! -=-
	 * As explained above, according to the FAQ, loading a collection of spectra "XSpectra"
	 * is slow, because at that moment Xcalibur attemps to compute all spectra whithin
	 * this is associated with significant overhead. (Specially in linux with Wine/Winelib)
	 *
	 * As an altenative, we will attemps to load *ALL* spectra at once
	 * and pray that the (32bit) OLE/COM system will be able to sustain the load.
	 *
	 * TODO get the rest of xenobol to compile as a Linux 64bits application
	 * TODO or at least find the way to run Xcalibur in a separate process as an external server
	 */
	CDispPtr spectcol;
		
	cout	<< endl
		<< "Attempting to load all spectra at once" << endl
		<< "Please be patient this may take a while" << endl
		<< endl << "loading ...";
		
	try {
		/**
		 * 1 single call to minimize the calls-to-OLE/COM overhead in Wine.
		 */
		int filnum = 0;
		WITH1(spectcol, detect, L"Spectra %d", filnum) { // <- this takes geological time
			cout << "done";
			try {
				int numspect = 0;
				dhCheck( dhGetValue(L"%d", &numspect, spectcol, L"Count") );
				cout << "\t" << numspect << " spectra" << endl;
	
				for (int s = 1; s <= numspect; s++) {
					try {
						unsigned short	realidx = 0, // plain simple zero-based index as in C (instead of 1-based as in Basic)
								scannum = 0; // absolute number of the scan (important when searchback a specific number)
						double retentiontime = -1;
						bool docentroid = false;
	
						cout << "\tspectrum " << s;
	
						/*
						 * Get information from collection Index (Scan number, retention time, etc.)
						 */
						if (SUCCEEDED( dhCallMethod(spectcol, L"IndexToScanNumber %d,%&hd,%&hd", s, &realidx, &scannum) )) {
							cout << " " << scannum;
						} else {
							cerr << "Warning, can't get an absolute scan number !!!" << endl;
						}
						
						if (SUCCEEDED( dhCallMethod(spectcol, L"IndexToRetentionTime %d,%&hd,%&le", s, &realidx, &retentiontime) )) {
							// Xcalibur produces minutes, we use seconds
							retentiontime *= 60.;
							cout << " at " << setprecision(4) << retentiontime << "s";
						} else {
							cerr << "Warning, can't get a retentiontime !!!" << endl;
						}
						
						/*
						* Get filter
						*/
						FilterParse* curfilt = NULL;
						
						try {
							WITH(filtercol, detect, L"Filters") {
								WITH1(filter, filtercol, L"ScanNumber %d", scannum) {
									char *str = NULL;
									dhCheck( dhGetValue(L"%s", &str, filter, L"Text") );
									cout << "\tfilter " << str;
									if (filters[0] != NULL)
										delete curfilt;
									filters[0] = curfilt = new FilterParse (str);
									cout << endl;
									SAFE_FREE_STRING(str);
								} END_WITH_THROW(filter);
							} END_WITH_THROW(filter);
						} catch (string errstr) {
							cerr	<< "Warning, can't get filter." << endl 
								<< errstr << endl
								<< "Trying to continue anyway..." << endl;
							continue;
						}
	
						cout << "\tType " << curfilt->ScanType() << " MS" << curfilt->MSLevel();
					
						/*
						* Process the data according the type...
						*/	
						if (curfilt->ScanType() == "") {
							// ..if we have a type
							cerr << "No Scan Type could be parsed from string. Skipping !!!" << endl;
							continue;
						} else if (curfilt->ScanType() == "SRM" || curfilt->ScanType() == "CRM") {
							// TODO build precursor scans out of Reaction monitor
							cout << "\tunsupported, skipping" << endl;
							continue;
						} else if (curfilt->MSLevel() == 2) {
							if (! dropms2) {
								WITH1(spect, spectcol, L"Item %d", s) {
									do_single_spectrum_ms2 (
										d, 0, 
										spect, scannum, retentiontime, 
										skipms2, docentroid
									);
								} END_WITH_THROW(spect);
							} else
								cout << "\tuser skipping all" << endl;
							continue;
						} else if (curfilt->MSLevel() == 1) {
							if (! dropms1) {
								WITH1(spect, spectcol, L"Item %d", s) {
									do_single_spectrum_ms1 (
										d, 0, 
										spect, scannum, retentiontime, 
										skipms1, docentroid
									);
								} END_WITH_THROW(spect);
							} else
								cout << "\tuser skipping all" << endl;
							continue;
						}
	
						// TODO implement other scans
						cout << "\tunknown, skipping" << endl;
						cout << endl;
	
					} catch (string errstr) {
						cerr	<< "Error getting spectrum collection" << endl 
							<< errstr << endl
							<< "Trying to continue with next spectrum..." << endl;
						continue;
					}
				}
				cout << endl;
			} catch (string errstr) {
				cerr	<< "Error while processing spectra" << endl
					<< errstr << endl
					<< "Aborting" << endl;
				return;
			}
		} END_WITH_THROW(spectcol);
	} catch(string err) {
		cout << "failed !" << endl;
		cerr	<< "System was unable to load a collection with all spectra inside." << endl
			<< err << endl;
		return;
	}
}




/**
 * This step decodes 1 single fragment spectrum
 * Normally, it is called by a loop in dospectra_ms2
 * BUT it can be used in other situations too (using special-case loops)
 *
 * General parameters :
 * @param d	detector we are currently to processing 
 * @param f	filter we are currently to processing 
 *
 * Parameters for spectrum (mostly data which is obtained outside the spectrum, from the collection) :
 * @param spect	OLE/COM dispatch object holding the spectrum to be analysed
 * @param scannum	unique number of the scan
 * @param retentiontime	
 *
 * Parameters how to process the data
 * @param docentroid	FIXME unimplemented - seems to be lacking in object internally of "XSpectrumRead" coclass.
 * @param skip	do we process the data ? or do we skip it
 */
void arthur::do_single_spectrum_ms2(const int d, const int f, CDispPtr spect, unsigned short scannum, double retentiontime, const bool skip, bool docentroid) {
	FilterParse* curfilt = (f>=0) ? filters[f] : NULL;
	
	/*
	* Get informations from header
	*/
	bool Centroided = false;
	try {
		WITH(head, spect, L"Header") {
			XSpectrumTypes type;
			dhGetValue(L"%d", &type, head, L"Type");
			if (type == XSpectrumCentroid)
				Centroided = true;
		} END_WITH_THROW(head);
	} catch (string err) {
		cerr << "Warning, can't get header - details" << endl
			<< err << endl;
	}

	if (Centroided) {
		cout << "\tcentroided";
	} else {
		if (Centroiding) {
			cout << "\tcentroiding";
			docentroid = true;
		} else {
			cout << "\traw";
		}
	}

	/*
	* Get information about parents
	*/
	int charge = -1;
	double mass = 0.;
	try {
		double m2 = 0.;
		int numparents = 0;
		unsigned parentscannum = 0;
		WITH(parentcol, spect, L"ParentScans") {
			dhCheck( dhGetValue(L"%d", &numparents, parentcol, L"Count") );
			cout << "\t" << numparents << " parents"<< endl;
			
			// HACK currently we only use parent n° 1 
			// TODO once we get support for MS3 and beyond in CaptainBol, we should upgrade this code
			if (numparents) WITH1(parent, parentcol, L"Item %d", 1) {
				dhCheck( dhGetValue(L"%d", &charge, parent, L"ChargeState") );
				cout << "\t\t\tcharge " << charge;

				dhCheck( dhGetValue(L"%e", &mass, parent, L"IsolationMass") );
				cout << "\tmass " << mass;

				dhCheck( dhGetValue(L"%e", &m2, parent, L"MonoIsoMass") );
				cout << " mono " << m2;

				dhCheck( dhGetValue(L"%d", &parentscannum, parent, L"ScanNumber") );
				cout << "\t#" << parentscannum;
				
				// memorise the association between current scan and its parent, used to make references in bol::runlcmsms later
				if (scannum < USHRT_MAX)
					parentmap[scannum] = parentscannum;

				// TODO we can access the parent's TIC using parent.Scan.Header.Tic
				cout <<  endl;
			} END_WITH_THROW(parent);
		} END_WITH_THROW(parentcol);
	} catch(string err) {
	//	cout << "\tno parents";
		cerr << "Warning, can't get parent informations !!!" << endl
				<< err << endl;
	}
	
	///
	/// At THIS point, we have a lot of information gathered. We can start building a spectrum
	///
	// *** EXPORT CAPTAINBOL
	
	// No parent ? Use info from Filter
	if (mass == 0.&& curfilt)
		mass = curfilt->Parent();
	// Charge mapper
	ChargeMapper chm;
	if (charge >= 1)
		chm.addCharge(charge);

	// spectrum !
	SpectrumLcmsms *sp = new SpectrumLcmsms(mass, chm.charges());
	if (sp == NULL) {
		cerr << "Critical Out of memory !!!" << endl;
		exit(255);
	}
	// other properties
	if (retentiontime >= 0.)
		sp->setRetentionTime(retentiontime);
	REFID r = { 0LL };
	r.th.scanid = scannum;
	r.th.detectid = d-1;
	sp->setRefid(r.refid);
	if (curfilt) {
		sp->setDescription(curfilt->Text());
		if (curfilt->Energy() > 0)
			sp->setCollisionEnergy(curfilt->Energy());
	}

	/// Store the spectrum
	run->msmsCollection().push_back(sp);// store into bol::runlcmsms
	parentsp2[scannum] = sp;	// memorize which scannum is at which address, used to make references in bol::runlcmsms later
	if ((! skip) || extractmap[scannum]) 
		getdata(sp, spect);
	else
		cout << "\t user skipping peaks" << endl;
}


/**
 * This step decodes 1 single parent spectrum
 * Normally, it is called by a loop in dospectra_ms1
 * BUT it can be used in other situations too (using special-case loops)
 *
 * General parameters :
 * @param d	detector we are currently to processing 
 * @param f	filter we are currently to processing 
 *
 * Parameters for spectrum (mostly data which is obtained outside the spectrum, from the collection) :
 * @param spect	OLE/COM dispatch object holding the spectrum to be analysed
 * @param scannum	unique number of the scan
 * @param retentiontime	
 *
 * Parameters how to process the data
 * @param docentroid	FIXME unimplemented - seems to be lacking in object internally of "XSpectrumRead" coclass.
 * @param skip	do we process the data ? or do we skip it
 */
void arthur::do_single_spectrum_ms1(const int d, const int f, CDispPtr spect, const unsigned short scannum, const double retentiontime, const bool skip, bool docentroid) {
	FilterParse* curfilt = (f>=0) ? filters[f] : NULL;

	/*
	 * Get informations from header
	 */
	bool Centroided = false;
	try {
		WITH(head, spect, L"Header") {
			XSpectrumTypes type;
			dhGetValue(L"%d", &type, head, L"Type");
			if (type == XSpectrumCentroid)
				Centroided = true;
		} END_WITH_THROW(head);
	} catch (string err) {
		cerr << "Warning, can't get header - details" << endl
			<< err << endl;
	}

	if (Centroided) {
		cout << "\tcentroided";
	} else {
/*		if (Centroiding) {
			cout << "\tcentroiding";
			docentroid = true;
		} else {
*/			cout << "\traw";
//		}
	}

	///
	/// At THIS point, we have a lot of information gathered. We can start building a spectrum
	///
	// *** EXPORT CAPTAINBOL

	// spectrum !
	SpectrumLcms *sp = new SpectrumLcms(retentiontime);
	if (sp == NULL) {
		cerr << "Critical Out of memory !!!" << endl;
		exit(255);
	}
	// other properties
	REFID r = { 0LL };
	r.th.scanid = scannum;
	r.th.detectid = d-1;
	sp->setRefid(r.refid);
	if (curfilt) {
		sp->setDescription(curfilt->Text());
	}

	/// Store the spectrum
	run->msCollection().push_back(sp);	// store into bol::runlcmsms
	parentsp1[scannum] = sp;	// memorize which scannum is at which address, used to make references in bol::runlcmsms later
	if ((! skip) || extractmap[scannum]) 
		getdata(sp, spect);
	else
		cout << "\t user skipping peaks" << endl;
}

/**
 * This function extract the peak data it self.
 * Note : This function throws exception
 *
 * @param sp bol::Spectrum into which to store data
 * @param spect XSpectrum object from which to get data
 */
void arthur::getdata(Spectrum< PeakIntensity* > *sp, CDispPtr spect) throw(string) {
	/**
	 * -=- !!!XDK!!! -=-
	 * Access to the data is a little bit complicated because
	 * the data is stored into a bi-dimensionnal array
	 *
	 * The array has several column
	 * - The first is the mass
	 * - The second is the intensity
	 * (the other aren't documented in the XDK)
	 *
	 * The array is handled using OLE's SafeArray.
	 * The DistHelper driver can only transmit them as Variant.
	 */
	
	// GET THE DATA : this is going to be tricky !!!
	
	// ...the function return a variant (and that's as far as DispHelper supports)
	VARIANT v;
	
	// get the variant !
	dhCheck( dhGetValue(L"%v", &v, spect, L"Data") );
	
	// ...we check if that variant holds an array of doubles (R8)
	if (V_VT(&v) != (VT_ARRAY|VT_R8)) {
		string err ("Bad data type : ");
		if (V_VT(&v) & VT_ARRAY) 
			err += "array of ";
		if (V_VT(&v) & VT_VECTOR)
			err += "vector of ";
		if (V_VT(&v) & VT_BYREF) 
			err += "pointer to ";
		err += "type " + boost::lexical_cast<string>(V_VT(&v) & VT_TYPEMASK);
		throw err;
	}
	
	// Got an array !
	SAFEARRAY *a = V_ARRAY(&v);
	// ...we need 2 Dimensions at least.
	if (SafeArrayGetDim(a) < 2) {
		string err ("Bad data array dimensions : ");
		err += boost::lexical_cast<string>(SafeArrayGetDim(a));
		err += "D";
		throw err;
	}
	// ...get the size in the 2 first dimensions
	LONG minc=-1,minp=-1,maxc=-1,maxp=-1,numc,nump; // min & number of peaks and data column per peaks.
	SafeArrayGetLBound(a,1,&minc);
	SafeArrayGetLBound(a,2,&minp);
	SafeArrayGetUBound(a,1,&maxc);
	SafeArrayGetUBound(a,2,&maxp);
	numc=1+maxc-minc;
	nump=1+maxp-minp;
	if ((minc<0)||(minp<0)||(maxc<0)||(maxp<0)||
		(nump <= 0)||(numc <= 0)) {
		string err ("Bad data array size : ");
		err += boost::lexical_cast<string>(numc);
		err += "x";
		err += boost::lexical_cast<string>(nump);
		err += "\t<";
		err += boost::lexical_cast<string>(minc);
		err += ":";
		err += boost::lexical_cast<string>(minp);
		err += ">-<";
		err += boost::lexical_cast<string>(maxc);
		err += ":";
		err += boost::lexical_cast<string>(maxp);
		err += ">";
		throw err;
	}
	cout << "\t" << nump << " peaks" << endl;

#ifndef _MSC_VER
	// Get a decent pointer NOW !
	double (*d)[numc];
	if (! SUCCEEDED( SafeArrayAccessData(a,reinterpret_cast<void**>(&d)) )) {
		string err ("Can't get pointer to data");
		throw err;
	}

	// Phew ! Now we're ready to work with the data
	if(sp) {
		// *** EXPORT CAPTAINBOL
		
		sp->resize(nump, NULL);
		for (int i = 0; i < nump; i++) {
			PeakIntensity* pk = new PeakIntensity(d[i][0],0,d[i][1]);
			if (NULL ==  pk) {
				cerr << "Critical Out of memory !!!" << endl;
				exit(255);
			}
			(*sp)[i] = pk;
		}
	} else {
		// dump to the screen
		for (int j = 0; j < nump; j++) {
			cout << "\t\t";
			for (int i = 0; i < numc; i++) {
				cout << "\t" << d[j][i];
			}
			cout << endl;
		}
	}
#else
	// According to microsoft :
	// C99 Variable length arrays are not currently supported in Visual C++. 

	// So wi will use a plain pointer to doubles
	double*d;
	if (! SUCCEEDED( SafeArrayAccessData(a,reinterpret_cast<void**>(&d)) )) {
		string err ("Can't get pointer to data");
		throw err;
	}

	// and do the ugly match pointer arithmetic ourselves.
	if(sp) {
		// *** EXPORT CAPTAINBOL
		
		sp->resize(nump, NULL);
		for (int i = 0; i < nump; i++) {
			PeakIntensity* pk = new PeakIntensity(d[i*numc+0],0,d[i*numc+1]);
			if (NULL ==  pk) {
				cerr << "Critical Out of memory !!!" << endl;
				exit(255);
			}
			(*sp)[i] = pk;
		}
	} else {
		// dump to the screen
		for (int j = 0; j < nump; j++) {
			cout << "\t\t";
			for (int i = 0; i < numc; i++) {
				cout << "\t" << d[j*numc+i];
			}
			cout << endl;
		}
	}
#endif
	SafeArrayUnaccessData(a);
	
//	if (! SUCCEEDED( SafeArrayDestroy(a) )) {
	if (! SUCCEEDED( VariantClear(&v) )) { // <- will recursively clear the Variant and all of its content
		string err ("Can't free memory ?!?");
		throw err;
	}
}


/**
 * Use the vectors "parent" & Co to populate references inside bol::runlcmsms
 *
 * parentmap - remembers, for a specific ms^2 spectra 'scannum0, which is the 'scannum' of the parent
 * parentms1 - a table which stores pointer to bol ms spectrum for a given XCalibur 'scannum'
 * parentms2 - a table which stores pointer to bol ms^2 spectrum for a given XCalibur 'scannum'
 */
void arthur::buildptr() {
	// *** EXPORT CAPTAINBOL
	for (SpectrumLcmsmsCollection::const_iterator i = run->msmsCollection().begin();
		    i != run->msmsCollection().end(); i++) {

		SpectrumPeakRefDualPtr<PeakIntensity>* rptr = NULL;
		double it = 0.;

		assert((*i) != NULL);

		// find back the scannum using refid
		REFID r = { (*i)->refid() };
		unsigned short scannum = r.th.scanid;
		unsigned short parentscannum = parentmap[scannum];
		Spectrum<PeakIntensity*>* sptr = static_cast < Spectrum< PeakIntensity* >* > (parentsp1[parentscannum]);
	
		// if the parent happens to be in the base....
		if (sptr != NULL) {
			cout << "scan " << scannum << " has parent " << parentscannum << " at address "
					<< static_cast<void*> (sptr) << endl;

			PeakIntensity* pptr = NULL;
			

			// find peak too (by mass)
			rptr = FindParentPeak(static_cast< SpectrumLcms* >(sptr), (*i), (*i)->moz(), &it);
		} else {
			cout << "scan " << scannum << " has no parent";
			
			// no parent : scan all precursor spectra
			rptr = FindParentSpectrum(run->msCollection(), (*i), (*i)->moz(), &it);
			if (rptr != NULL) {
				REFID r = { rptr->spectrum()->refid() };
				parentscannum =  r.th.scanid;
				
				cout << "\tdiscovered " << parentscannum << " at addess "
					<< static_cast<const void*> (rptr->spectrum());
			}
			cout << endl;
		}
		
		// mutate the intensity to the new value
		if (it > 0.)
			(*i)->setIntensity(it);

		if (rptr != NULL) {
			// we have a reference to a MS^1 spectrum & peak....
			SpectrumPeakRefPtr baseptr(rptr);
			(*i)->setSpectrumPeakRef(baseptr);
	
			// if everything is okay, now we can travel from the current spectrum to the parent intensity
			// (and both should be the same)
			if (rptr->peak()) // check if there's one
				assert((*i)->spectrumPeakRef()->peak()->intensity() == it);
		}
	}
}


/**
 * This function builds a vector of boolean specifying which 'scannum' have to be extracted
 *
 * @param detector current detector
 * @param ptr pointer to a vector to 64bits refid
 */
void arthur::initextract_map(unsigned short d, ExtractList *ptr) {
	extractmap.clear(); extractmap.resize(USHRT_MAX,false);
	
	if (ptr)
		for (ExtractList::iterator i = ptr->begin(); i < ptr->end(); i++) {
			// use the union to convert the 64bits refid into numbers we need
			REFID rf = { *i };
			if (rf.th.detectid == d-1)
				extractmap[rf.th.scanid] = true;
		}
}


}
}
